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Seismic images of Earth's interior have revealed two continent-sized anomalies 

with low seismic velocities, known asthe large low-velocity provinces (LLVPs), 

in the lowermost mantle’. The LLVPs are often interpreted as intrinsically dense 
heterogeneities that are compositionally distinct from the surrounding mantle?. 

Here we show that LLVPs may represent buried relics of Theia mantle material (TMM) 
that was preserved in proto-Earth’s mantle after the Moon-forming giant impact’. Our 
canonical giant-impact simulations show thata fraction of Theia's mantle could have 
been delivered to proto-Earth's solid lower mantle. We find that TMM is intrinsically 
2.0-3.596 denser than proto-Earth's mantle based on models of Theia's mantle and the 


observed higher FeO content of the Moon. Our mantle convection models show that 
dense TMM blobs with a size of tens of kilometres after the impact can later sink and 
accumulate into LLVP-like thermochemical piles atop Earth’s core and survive to the 
present day. The LLVPs may, thus, bea natural consequence of the Moon-forming 
giant impact. Because giant impacts are common at the end stages of planet accretion, 
similar mantle heterogeneities caused by impacts may also exist in the interiors of 
other planetary bodies. 


Terrestrial planet formation concludes with a number of giant impacts 
among Moon- and Mars-sized planetary embryos’. The well-studied 
giant-impact scenario involves a protoplanet, Theia, colliding with 
proto-Earth. Earth’s Moon subsequently formed from orbitally bound 
impact debris*. Although it is challenging to reconcile the observed 
lunar isotopes and volatiles®’ with the giant-impact hypothesis?, 
this Moon-formation hypothesis accounts for several key features of 
the Earth-Moon system, including the current angular momentum, the 
Moon’s small core and the high mass of the Moon compared with the 
Earth?°. However, direct evidence for the existence of Theia remains 
elusive. Here, we demonstrate that a fraction of Theia mantle mate- 
rial (TMM) could have entered the largely solid lower layer of the 
post-impact Earth’s mantle during the canonical Moon-forming impact. 
This mostly molten TMM could later have solidified and sunk to Earth’s 
lowermost mantle and now constitute the seismically observed large 
low-velocity provinces (LLVPs)! in the present day (Fig. 1). 

Previous hypotheses for the origin of the LLVPs, such as remnants 
of Earth's early differentiation?? and accumulations of subducted 
oceanic crust’, specifically involve Earth's internal processes after 
the Moon-forming giant impact. However, lines of evidence suggest 
that the LLVPs preserve primordial volatiles that predate the giant 
impact?! and these volatiles cannot easily be explained by Earth's 
ensuing differentiation after the Moon's formation. Since Theia pre- 
dates the Moon, its mantle may havebeen ableto retain nebular gases" 
that were captured due to its large mass before the dissipation of its 


protoplanetary disc”. Thus, a TMM-origin of the LLVPs may explain 
the presence of nebular components in some ocean-island basalts 
(OIBs)! that are caused by mantle plumes that entrain materials from 
the LLVPs“. This mechanism can also explain the similarity of noble-gas 
isotopic ratios between OIBs and lunar mare basalts». However, sev- 
eral additional physical and chemical criteria must be met for TMM 
to be the origin of LLVPs. In what follows, we combine evidence from 
hydrodynamic simulations of the Moon-forming giant impact, ther- 
mal evolution modelling and mantle convection simulations to show 
the following. (1) The giant impact could have produced a two-layer 
structure in the Earth's mantle, with an upper molten layer and a lower 
solid layer. (2) The solid lower layer could have directly incorporated 
some mostly molten TMM. (3) This molten TMM could solidify, sink and 
accumulate atthe Earth's core- mantle boundary (CMB). (4) Accumu- 
lations of TMM on the CMB could survive more than 4.5 billion years 
without being advected or eroded away from the CMB. (5) The volume 
andseismic properties ofthe TMM accumulations are consistent with 
those ofthe LLVPs. 

We investigated the dynamics and thermal state of Earth's mantle 
and the TMM during the giant impact using two different hydrody- 
namical simulation methods with the latest equations of state and 
at unprecedentedly high resolution". Using a meshless finite mass 
(MFM) method’, recent work? has shown that a Moon-forming giant 
impact could lead to a stratified Earth mantle characterized by a com- 
positional change and entropy jump near a depth of approximately 
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Fig.1|Schematic diagram illustrating the giant-impact origin of the LLVPs. 
a, Pre-impact: Theia collides with the proto-Earth. b, Post-impact: a stratified 
mantle forms in the Earth. There is a boundary at a depth of approximately 
1,400 km within the present-day mantle (marked by the dashed white line). 
The mantle of post-impact Earth (blue) below the boundary is largely solid, 

in contrast with its molten upper layer. Most of the TMM (red) is incorporated 
in this upper layer, but a portion of the TMM enters the lower layer. Data and 


1,300 km (referenced to Earth’s present-day surface radius). From 
these results”, we found that a significant amount of TMM, 0.017- 
0.026 Earth masses (Mọ) (Extended Data Table 1), could have entered 
the lower mantle layer, which is comparable with mass estimates of 
approximately 0.01-0.06 M; for the LLVPs“”°. However, both thermo- 
dynamics and material mixing are sensitive to the setup of numerical 
models, including equations of state, numerical resolution and the 
hydrodynamical method (Methods). To further assess the delivery of 
TMM tothe deep mantle, we performed new MFM impact simulations 
with higher numerical resolution and updated equations of state! 
(Methods). The new results reproduce a similar layered structure in 
the post-impact Earth as in ref. 19, with the upper layer dominated 
by melt and vapour, and a rheologically solid lower layer (Extended 
Data Fig. land Supplementary Video 1). The total mass of TMM within 
the lower layer is approximately 0.01 Earth masses, and 15% of the 
lower-layer TMM can be solid (Extended Data Table 2). We also per- 
formed ultra-high-resolution smoothed particle hydrodynamics (SPH) 
simulations of giant impacts” to corroborate the results (Methods). 
The high-resolution SPH simulations also demonstrate, as with the 
MFM simulations, that post-impact Earth could have retained a solid 
lower layer with some solid and molten TMM material (Supplementary 
Video 2). We found that as the numerical resolution was increased, 
the amount of TMM mass in Earth's lower layer and its solid fraction 
(up to 30%) increased considerably (Extended Data Table 3). TMM may 
mix with Earth's mantle at a scale smaller than our model's resolution; 
therefore, the size and dispersal pattern of TMM in the lower layer is 
not yet well constrained by our giant-impact models. Although experi- 
ments show efficient fragmentation of a body of dense fluid sinking 
through less dense fluid, they also show that the fragments remain 
tightly clustered and sink quickly at rates determined by their collec- 
tive behaviour and not by their individual Stokes velocities”. Future 
laboratory experiments and higher-resolution simulations are needed 
to investigate the small-scale dynamics during the giant impact. 

The composition of the TMM determined its intrinsic density, and 
thus the nature of its long-term advection in Earth’s mantle. In most 
Moon-formation impact scenarios, the Moon is preferentially com- 
posed of Theia’s mantle’. Given the higher FeO content of the Moon's 
mantle (more than 10 wt%) compared with Earth’s mantle (less than 
8 wt%)”, Theia’s mantle may also have been iron-rich compared to 
Earth's. Some ofthe enrichment of the Moon's FeO may have resulted 
from post-impact mechanisms", not from the high FeO content of 
Theia. Nevertheless, under different impact scenarios, the FeO content 
of Theia's mantle is estimated to be 11-40 wt%, based on Si isotopic 
constraints”. Recently, this model was extended” to also explain the 
very low2H/H (deuterium/hydrogen) ratio (3.91 x 10?) of inferred lunar 
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stratification model from ref. 19. c, Present day: after approximately 4.5 Gyr, the 
TMM inthe lower layer sinks to the CMB andis shaped into two thermochemical 
piles that represent the LLVPs. The rest of the TMM inthe upper layer undergoes 
significant mixing with the Earth’s mantle, whichis compositionally similar to 
the protolunar disc that formed the Moon. The relative sizes of the Earth and 
Moon andthe distance between them are not to scale in c. 


interior reservoirs”. This model predicts that Theia’s mantle contained 
approximately 13-18 wt% FeO (ref. 12). 

Based on the results of Desch and Robinson”, we implemented three 
composition models for the TMM with a fixed molar ratio of Mg/Si=1 
and a variable FeO content of 13, 15 or 17 wt% (Extended Data Table 4) 
to estimate the TMM density. We calculated the phase equilibria and 
density of TMM along Earth's geotherms using the Perple_X””78 code 
and a thermodynamic dataset”. We then compared the results with 
the densities of bulk silicate Earth (BSE)”. The three TMM composition 
models and the BSE yielded similar mineral assemblages of olivine, 
garnet and pyroxene under the pressure and temperature (P-7) condi- 
tions of the Earth’s upper mantle, which transform into bridgmanite 
(Mg,Fe,Al)(Si,Al)O,, magnesiowüstite (Mg,Fe)O and davemaoite CaSiO, 
under Earth's lower-mantle P-7 conditions (Extended Data Figs. 2 and 
3). We found that the TMM is approximately 2.0-3.5% denser than the 
BSE (Fig. 2and Extended Data Fig. 4) throughout Earth’s lower mantle. 

After the impact, mantle convection is expected to develop first from 
the outer part of Earth’s thermally stratified mantle??, which may then 
have propagated downwards and mixed with the lower-mantle TMM. 
However, this viscous process involves a prolonged time frame during 
which the dense TMM inthe lower layer may already have sunk deeper 
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Fig. 2 | Density profiles of the TMM and the BSE as a function of pressure. 
The density was calculated using Perple X? with thermodynamic data from 
ref. 29. The upper and lower bounds of the density profiles are, respectively, 
basedongeotherms from Brown and Shankland” and Stacey. 
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Fig.3| Theformation of LLVP-like thermochemical piles from intrinsically 
dense TMM. a, Initial composition field ina geodynamic model in which the 
TMM is introduced into the model domain from a depth of 1,400 km to the 
CMBasrandomly distributed spheres. b, Temperature field of the model at 
4.5 Gyr. c, Composition field at 4.5 Gyr showing the formation of LLVP-like 
thermochemical piles dueto the accumulation of TMM. The composition fields 
in panels a and c show pure background mantle material (0), pure TMM (1) ora 


to the CMB (Extended Data Fig. 5), thus avoiding significant stirring 
and mixing with materials in the upper layer. By the time it reached the 
CMB, the molten fraction of the TMM may also have mostly solidified 
due to cooling from the surroundings plus adiabatic compression. Our 
scaling analysis shows that a molten TMM blob would cool rapidly until 
it experiences a rheological transition. We estimated that the solidifi- 
cation timescale of a 50 km TMM blob can be as short as 1,000 years 
(Methods). Smaller blobs of TMM would solidify more quickly, and it 
is intuitive that they may sink less efficiently. However, their sinking 
tendency is heavily influenced by whether they are evenly dispersed 
within the background mantle matrix or remain clustered together” 
as well as the vigour of the convective flow in the background. To pre- 
cisely determine the initial state, size, distribution and differentiation 
process of lower-layer TMM right before the long-term convection of 
Earth’s solid mantle would require a combination of different equations 
of state for Theia and proto-Earth, unparalleled resolution for impact 
simulations and multiple-phase flow dynamics, which all merit future 
investigation. Inthe following, we assume a relatively short solidifica- 
tion timescale and treat the initially mixed molten and solid TMM in 


mixture (intermediate values). d, Seismic shear velocity perturbations (6V,) in 
the geodynamic model at the lower mantle depths. 5V, was calculated for each 
point by subtracting the average horizontal shear velocity and then dividing by 
this average. An adiabatic temperature gradient of 0.35 K km was added to the 
temperature field to calculate the shear-wave velocity. The density and V, were 
calculated using Perple_X?”8 with the thermodynamic dataset from ref. 29. 


the lower mantle as all rheologically solid to explore their long-term 
dynamics using numerical mantle convection simulations. 

Previous geodynamic modelling experiments have shown that chemi- 
cal heterogeneities with an intrinsic density a few percent higher than 
the surrounding mantle can sink to and accumulate at the CMB”! * and 
that these heterogeneities can later be shaped into spatially isolated 
thermochemical piles that resemble the seismically observed LLVPs**. 
We designed geodynamic simulations to further investigate the advec- 
tion of intrinsically dense TMM in Earth’s mantle throughout Earth’s 
history (Methods). In our reference geodynamic model, the TMM is 
2.5% intrinsically denser than the surrounding mantle and is initially 
introduced to the lower half of the model domain as randomly distrib- 
uted spheres with a fiducial radius of 50 km and a cumulative volume of 
approximately 4% of the two-dimensional (2D) model domain (Fig. 3a). 
The model starts with a uniformly hot temperature” to simulate the 
early hot mantle after the impact, but our results are not sensitive to 
this choice, as discussed later. After this initial state, the TMM quickly 
descends to the CMB andis later shaped by mantle convection into two 
spatially isolated thermochemical piles after 4.5 Gyr (Supplementary 
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Video 3). Using the temperature field and the composition field from 
this model (Fig. 3b,c), we calculated the seismic shear-wave velocity of 
the TMM piles and found that it is approximately 1-5% lower than the 
surrounding mantle (Fig. 3d), in agreement with the observed seismic 
reduction of the LLVPs/?955?? In the model, the volume of the piles is 
approximately 3.4?6 of the mantle, which aligns with the LLVPs cur- 
rently occupying a comparable volume percentage (approximately 
4%) of Earth's mantle*°. 

We additionally tested eight models with different densities 
(buoyancies) for the TMM (1.25%, 3.5% or 5% denser), a smaller radius 
for TMM blobs (25 km), radioactive heating of the TMM blobs, higher 
temperature dependence of the viscosity (x10? viscosity range owing 
to changes intemperature), periodic side-boundary conditions anda 
relatively cold initial temperature (Supplementary Videos 4-11). The 
physical parameters anda full list of geodynamic models are provided 
in Extended Data Tables 5 and 6, respectively. All of these models are 
characterized by similar LLVP-like piles above the CMB after 4.5 Gyr, 
except for one case in which the TMM was only 1.25% denser than the 
surrounding mantle and had a radius of 50 km and another case in 
which the TMM was 2.5% denser compared with the background mantle 
and had a small radius of 25 km. In all other cases, the higher density 
and larger initial size of TMM blobs increased the possibility that the 
TMM would sink and accumulate on the CMB and coherently persist 
after 4.5 Gyr. 

Giant accretionary collisions can create widespread melting in the 
form of partial or whole silicate magma oceans *?. Homogenization 
may take place during or after an impact. During solidification of the 
magma ocean, gravitational instability and rapid mixing can occur 
within the solidified region, due to convection triggered by a supera- 
diabaticthermal profile”. Alternatively, gravitational instability can 
be triggered by the Fe enrichment of the upper layers during crystal- 
lization*^*, Recent studies have also shown that convection within the 
solidified mantle could become more efficient as the magma ocean 
crystallizes’. However, the TMMs we focus on here are incorporated 
into the solid lower layer, which has a lower potential temperature 
than the overlying molten layer. Thus, the solid and relatively cold 
lower layer would prohibit whole-mantle convection shortly after the 
impact”, which limits mass exchange with the upper layer. Being denser 
than the solid background mantle, the TMM blobs considered here 
are likely to sink in the lower layer*®. Similarly, the dense TMM in the 
upper molten layer of post-impact Earth, ifit has not been completely 
mixed into the proto-Earth???, may sink as well. A vigorously convecting 
magma ocean in this upper layer may lead to a large extent of mixing 
with the solidified mantle”, but the effects of planetary rotation could 
promote the development of even large-scale mantle heterogenei- 
ties during magma ocean crystallization”. Thus, the approximately 
0.014-0.026 M, mass of TMM in the lower layer constrained in ref. 19 
and this work isa conservative estimation of the total TMM entering the 
lower mantle. If the TMM remains partially molten, it may eventually 
accumulate at the CMB and constitute a dense basal magma ocean. A 
dense melt layer was assumed previously. It has been well established 
that it would fractionally crystallize as stable solids”, which are later 
shaped by mantle convection into LLVP-like thermochemical piles. 

The preservation of large-scale TMM structures in Earth’s lowermost 
mantle hasimportant implications for the composition and evolution of 
Earth’s interior, the origin of the Moon and the early accretion history of 
the Solar System. The presence of preserved TMM piles in Earth's lower 
mantle implies that Earth's lowermost mantle contained large-scale het- 
erogeneities from the beginning of Earth’s evolution. The subsequent 
approximately 4.5 billion years of mantle convection did not erase such 
structures formed during Earth’s primary accretionary period, con- 
sistent with evidence from the noble gases in plume-derived OIBs!??, 
Our model, hence, suggests that constraining the composition and 
evolution of Earth's mantle requires consideration of heterogeneous 
accretions during the stage of giant impacts. Furthermore, the survival 
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of dense basal TMM accumulations may suggest that any pre-existing 
ancientreservoirs (before the Moon-forming impact)? may also have 
been preserved, which could explain the observed multiple isotopic 
signatures that require pre-impact reservoirs”. 

Present-day LLVPs may be a combination of TMM and other com- 
positional heterogeneities, for example, former subducted oceanic 
crust*^, However, the TMM may not have fully mixed with other com- 
ponents, which aligns with isotopic evidence from some OIBs (having 
both high ?He/*He and anomalous ?W) and suggests the preservation 
of some ancient, deep-mantle, primordial reservoirs that were least 
modified by recycled crust. Although some of these signatures may 
be attributed to core-mantle interactions in the Earth???" che fact that 
similar signatures of primordial noble gases are found in KREEP-bearing 
mare basalts” suggests the signatures in both bodies may derive from 
TMM. Further studies about the chemical and isotopic compositions 
of LLVP-related plume samples, and comparisons to lunar mare basalts 
may be used to further test the hypothesis. Finally, given how com- 
mon giant impacts were during the late stages of planetary accretion??, 
compositional heterogeneities from large impactors may be preserved 
in other solid planets as well. 
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Methods 


Giant-impact simulations 

Simulations of the giant impact are subject to uncertainties arising from 
the hydrodynamic methods? the treatment of contact discontinui- 
tiesandfreesurfaces””, theequationsof state’ and thenumerical reso- 
lution””**, The SPH method is widely used in giant-impact simulations’. 
However, standard SPH schemes and common forms of artificial 
viscosity, which are used to maintain differentiability across shocksin 
SPH, allow for the triggering of numerical viscosity that can suppress 
turbulence and mixing and cause overheating”. Methodological 
developments to limit the application of artificial viscosity are ongo- 
ing, including revisions to the SPH equations®, developing new forms 
of artificial viscosity (for example, ref. 69) or using methods without 
artificial viscosity altogether". Simply increasing the resolution also 
reduces its fractional effect on the whole simulation. Here we analyse 
models from ultra-high-resolution SPH simulations" and a new mesh- 
less method”. Despite the differences in the hydrodynamic methods, 
equations of state and resolution, we derived qualitatively similar 
results, namely, a solid lower mantle in the proto-Earth after a canoni- 
cal impact with approximately 0.01 M, of Theia material implanted 
into it. 


Giant-impact models with SPH simulations 

TheSPH simulations were run using the open-source hydrodynamics 
and gravity code SWIFT”, using 10? SPH particles (and lower-resolution 
comparisons). The initial conditions and full details are provided in 
ref. 17. The initial conditions were created using the publicly avail- 
able WoMA and SEAGEN packages”. Our numerical resolution is 
considerably greater than the 10°-10° particle resolution often used 
in standard Moon-formation studies, which enables the detailed 
study of the relatively small mass of TMM implanted into the deep 
Earth mantle. 

To model the core and mantle materials, we use the updated ANEOS 
Fe,<Si,, and forsterite equations of state. These allow direct exami- 
nation of the phase state of TMM and other material, in addition to 
checking the overall change in entropy. The physical parameters for 
the impact simulations are the same as in ref. 17. 

The simulations involve collisions between a low-mass (approxi- 
mately 0.1 Mọ) Theia and a large proto-Earth (0.9 Mẹ), atimpactangles 
of 45° and impact speeds at contact near the mutual escape speed, but 
with seven different model resolutions (Extended Data Table 3). All 
simulations produced a similar layered structure for the post-impact 
Earth with a mostly solid lower mantle. However, the mass and solid 
fraction of TMM implanted into the lower layer are not adequately 
resolved for numbers of particles below 107. Even with 10° particles, the 
precise values have not fully converged, although the overall behaviour 
is consistent. 

Although the simulations conducted in this study usedlow impactor 
masses of approximately 0.1 Me, we anticipate that the preservation of 
TMM in Earth's deep mantle may be also possible if Theia has a larger 
mass, as more TMMis likely to be implanted at greater initial depths and 
the FeO content (and thus density) of TMM could be higher”. There- 
fore, the Theia origin for the LLVPs that we consider here may be not 
restricted to the canonical-like impact scenarios, but the post-impact 
composition, thermal state of the lower mantle and mantle mixing in a 
wider range of impact conditions? would need to be investigated in 
detail, including with robust numerical treatments, to directly assess 
this claim. 


Giant-impact models in MFM simulations 

The MFM simulations were carried out with the open-source code 
GIZMO”. The method captures shocks with a Riemann solver and is, 
thus, free of artificial viscosity. MFM captures subsonic turbulence 
better than traditional SPH at a comparable resolution? and, thus, 


mixing in giant-impact simulations?. The MFM simulations cover 
alarger parameter space than the expensive ultra-high-resolution 
SPH simulations. The simulations described in ref. 19 used the 
ANEOS/M-ANEOS”* equations of state with iron comprising the core 
and forsterite comprising the mantle. It employed 5 x 10° particles. To 
explore the robustness of our results with respect to the equations of 
states and resolution, we then adopted the latest equations of state 
fromref.16 and ran new simulations, as listed in Extended Data Table 2, 
with 2 x 10° or 10 x 106 particles. 

Thanks to the new equations of state with explicit phase informa- 
tion, wecan identify the rheological transition in the mantle material 
(Extended Data Fig. 1). The radius of this rheological transition lies 
slightly beneath the entropy jump reported in ref. 19. Such entropy 
jumps are due to focused shock heating and concurrent supersonic 
flows in the upper mantle and are coincident with a compositional 
change”. More than 0.01 M, of TMMis consistently found in the lower 
Earth mantle, which is largely solid. 


Thermal evolution model of Theia mantle blobs 

We modelled the thermal evolution of the approximately 9096 of the 
TMMblobsinthelowerlayerthat were non-solid. These non-solid TMM 
blobs havea potential temperature of 5,000-8,000 K after the impact. 
The results of high-pressure melting experiments suggest that such 
blobs are entirely molten, even in the deep mantle”. They are, thus, 
likely to be cooled by thermal convection until the melt fraction 
becomes low enough to experience the rheological transition” ”*. We 
estimated the cooling timescale using a scaling analysis of convective 
heat flux. The heat flux out of the upper hemisphere F [W] is roughly 
given by 


where k; is the thermal conductivity, riy and Trymare the radius and 
temperature of Theia mantle blobs, 7,,is the temperature of the sur- 
rounding mantle, Rais the internal Rayleigh number and Ra, (approx- 
imately 1,000) is the critical Rayleigh number. Assuming a melt 
viscosity of  - 100 Pa s and a heat capacity of 1,000 J kg ! K*, the cool- 
ing timescale is predicted to be shorter than 1,000 years for a 
50-km-sized blob, regardless of the initial temperature. During this 
period, the TMM melt may also migrate through the surrounding man- 
tle by melt percolation. However, according to the giant-impact mod- 
els, the temperature ofthe deep lower mantle is below the solidus, and 
thus, achannelized melt flow would likely freeze soon after leaving the 
TMM blobs. Furthermore, even if the TMM melt were able to raise the 
temperature of the surrounding mantle to melt it, the porosity expected 
would be minimal, and therefore, percolation through such a mostly 
solidified mantle would be negligible. Additionally, given that the TMM 
blobs are embedded in the ductile lower mantle, crack formation via 
fracturing is not an effective mechanism for melt migration. 

The subsequent evolution of rheologically solid TMM would again be 
governed by percolation, which separates the melt and liquid phases 
within the blob. Depending on the density contrast between the melt 
and solid, the melt may migrate either upwards or downwards” *!, 
creating some degree of fractionation. Regardless of the direction, 
however, the melt would be expected to solidify as it reaches the rim 
ofthe Theia mantle blobs. The melt would also solidify by adiabatic 
compression as the blobs sink through the ambient mantle. Consid- 
eringthe estimated short percolation timescale in comparison to the 
long timescale of the viscous flow, it is plausible that the blobs would 
solidify before they reach the CMB. 


Setup of the mantle convection model 
Here, we provide the governing equations that are used to solve the 
mantle convection models in our modified 2D Citcom code. Complete 


descriptions are available from ref. 82 and references therein. We 
performed the geodynamic calculations by solving the following 
non-dimensional equations for conservation of mass, momentum 
and energy under the Boussinesq approximation: 


V-u=0 (1) 
-VP+V- (ne) - Ra(T- BC)z (2) 
oT +(uV)T=VT+Q (3) 


where uis the velocity, P isthe dynamic pressure, nis the viscosity, € is 
the strain rate, Ra is the Rayleigh number, 7 is the temperature, and 
BandCarethe buoyancy number and composition, respectively. zis the 
unit vector in the vertical direction, Cis time and Qis the internal heat- 
ing rate. a and x are the thermal expansivity and thermal diffusivity, 
respectively. 

The advection of the composition Cis given by: 


oC + (u-V)C-0 (4) 


The viscosity is defined as: 


n=N,N.explA(0.5- 7)] (5) 


where A is the activation coefficient for the temperature dependence 
of the viscosity. n isa parameter that determines the intrinsic viscosity. 
The viscosity prefactor 1, is 1.0 in the upper mantle and 50.0 in 
the lower mantle, resulting in a x50 viscosity increase across the 
660 km discontinuity. We used an activation coefficient of A = 9.21 in 
most cases, but a value of A = 11.51 was also explored to test an even 
higher temperature-dependent viscosity range owing to changes in 
temperature. 

The simulation was computed in a 2D Cartesian geometry with an 
aspect ratio of 5:1. To solve the conservation equations, we used a 
modified version of the convection code Citcom?, which includes 
thermochemical convection and both composition- and temperature- 
dependent rheology. As is appropriate for 2D Cartesian models®, we 
excluded internal heating in most calculations, but models with radio- 
active heating of the TMM were also tested. The initial temperature 
was 0.72 (non-dimensional) everywhere with small perturbations, and 
wealsotested a model with a different initial temperature. The top and 
bottom surfaces were isothermal with T= O and T- 1, respectively. All 
boundaries were free-slip. The side boundaries were reflective and 
insulating, but we also tested a model with periodic side boundaries. 
Thecompositional field was advected using approximately 7.9 million 
tracers with the ratio-tracer method™. The diffusion time was dimen- 
sionalized using t= kE, where t is dimensional time, t’ is non- 
dimensional time, his the thickness of the mantle (here, h = 2,890 km) 
and xis the thermal diffusivity (here, K = 1 x 10 5m? s). 


Data availability 

All data and parameters are available in the main text or the supple- 
mentary materials. The data that support the findings of this study 
are also available at https://doi.org/10.6084/m9.figshare.24013776.v1. 
Source data are provided with this paper. 


Code availability 


The author's modified 2D Citcom code used in this study is avail- 
able from https://figshare.com/projects/Yuan Li 2022 NG/129185. 
The GIZMO code is made available at http://www.tapir.caltech. 


edu/-phopkins/Site/GIZMO.html. SWIFT is publicly available at 
http://swiftsim.com. WoMais publicly available at https://github.com/ 
srbonilla/WoMa, or the Python module can be installed directly with 


pip (https://pypi.org/project/woma/). 
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Extended DataFig.1|Entropy profile (blue, inJkg "K ?) of mantle material 
in the post-impact Earth for our impact model using the meshless finite 
mass (MFM) method”. The red curve shows the liquidus of forsterite and 
the rheological transition of the mantle was marked by the orange curve where 
the melt fraction above that depth becomes larger than 4096*255, 
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Extended Data Fig. 2 | Phase diagrams ofthe bulk silicate Earth (a), Theia 1 
(b), Theia 2 (c) and Theia_3 (d) with geotherm from ref. 57. The FeO contents 
of Theiaare13 wt% (Theia 1), 15 wt% (Theia 2), and 17 wt% (Theia 3), respectively. 
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Phase equilibria were calculated using Perple X??? with thermodynamic data 
fromref.29.St:stishovite, Fp: ferropericlase, Ring: ringwoodite, Wad: wadsleyite, 
Ol: olivine, Cpx: clinopyroxene, Brg: bridgmanite, Gt: garnet, CaPv: davemaoite. 
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Extended Data Fig. 3 | Phase diagrams of the bulk silicate Earth (a), Theia 1 
(b), Theia 2 (c) and Theia_3 (d) with geotherm from ref. 58. The FeO contents 
of Theiaare13 wt% (Theia 1), 15 wt% (Theia 2), and 17 wt% (Theia 3), respectively. 
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Phase equilibria were calculated using Perple X??? with thermodynamic data 
from ref. 29. Fp: ferropericlase, Wad: wadsleyite, Ol: olivine, Cpx: clinopyroxene, 
Brg: bridgmanite, Gt: garnet, CaPv: davemaoite. 
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and the bulk silicate Earth as a function of pressure. Phase equilibria were 
calculated using Perple X?75 with thermodynamic data from ref. 29. 
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Extended Data Fig. 5| One numerical experiment showing that dense TMM blobs are placed in the lower mantle (b). After 27.36 Myr, the TMM blobs 
sinks tothe CMB before upper mantle materials mix with lower mantle reach the CMB (d) whereas there is little mixing between the upper mantle 
materials. a-d, Snapshots of the temperature fields (a, c) and compositional and lower mantle (c). 
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Extended Data Table 1| The depth of layer boundary and the enclosing mass of Theia mantle material below the boundary for 
the canonical Moon-forming giant impact simulations from ref. 19 


Impact model iin of layer Mass of Theia's mantle 
oundary (km) in lower layer (Mo) 

2 1,271 0.017 

3 1,222 0.026 

5 1,313 0.023 

6 1,313 0.024 
Average 1,280 0.023 
lo 43 0.0038 


Here the layer boundary is defined by a compositional change, i.e., the fraction of TMM in the mantle". We measure the fraction of enclosed mass below the layer boundary and then invert the 
enclosed mass radius relation in the PREM model? to get the depth of the layer boundary in present-day Earth. 


Extended Data Table 2 | Results from the MFM simulations with updated versions of the ANEOS equations of state'? 


Entr Mass of TMM solid Mass of TMM solid 
Run eR. TMM in fraction Rheological TMM in fraction 
(R H lower in lower boundary lower layer in lower 
® layer (Mẹ) layer (%) (Mo) layer (96) 
1 0.85 0.016 8.2 0.74 0.0088 15 
2 0.83 0.014 8.1 0.75 0.0084 5 
3 0.81 0.013 11. 0.73 0.0082 16 


The radius above which the mantle transit to a rheological liquid state (melt fraction >40%), the mass of TMM below that boundary, and the fractions of it that are solid. For comparison, we also 
listed the results using the boundary defined by entropy jump as in to ref. 19. These canonical impacts involve a low-mass (~0.12 M;) Theia and a large proto-Earth (0.9 M;) at -45? impact angles, 
which are similar to ref. 19 but with enhanced resolution and an updated equation of state. The third MFM simulation has a particle count that is five times higher compared to the other two cases. 
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Extended Data Table 3 | Results from the SPH simulations 


Run Number of Rheological Mass of TMM in TMM solid fraction 


particles boundary (Rẹ) lower layer (Mẹ) in lower layer (%) 
1 105 0.779 3.88x10? 0.0 
2 1055 0.804 1.52x104 0.0 
3 10° 0.811 4.62x104 0.6 
4 1065 0.805 6.57x10^ 0.3 
5 107 0.805 1.17x10? 19.6 
6 1075 0.814 1.54x10? 31.1 
7 108 0.823 2.59x10? 29.3 


The radius above which the fraction of solid mantle material has dropped to 60%, the mass of TMM below that boundary, and the fraction of it that is solid, for different simulation resolutions. 


Extended Data Table 4 | Major elemental compositions of the mantle of Theia and proto-Earth, and the bulk silicate Earth 
(BSE) (wt?6) used in our thermodynamic modeling 


Composition SiO» AbO; FeO MgO CaO Total 
Theia 1 47.72 4.03 13.00 32.01 3.24 100.00 
Theia 2 46.52 4.03 15.00 31.21 3.24 100.00 
Theia 3 45.32 4.03 17.00 30.41 3.24 100.00 

Proto-Earth 57.25 4.03 5.00 30.48 3.24 100.00 


Bulk silicate Earth 47.95 4.03 7.79 36.98 3.24 100.00 


The BSE composition is normalized to 100% from ref. 23. 
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Extended Data Table 5 | Physical parameters for mantle convection models 


Parameter Symbol Value 
Input Rayleigh number Ra 1x108 
Earth radius R 6,371 km 
Core—mantle boundary radius r 3,503 km 
Mantle thickness h 2,868 km 
Thermal expansivity Qo 1x10°/K 
Thermal diffusivity Ko 1x10 m?/s 
Gravitational acceleration 9o 9.8 m/s? 
cape Dicen nien ar 250K 
Background mantle dp 3,300 kg/m? 


reference density 


Reference viscosity No 5.2 x10?! Pas 


Extended Data Table 6 | Full list of performed mantle convection models 


Buoyancy Side 
Gas Ra number Activation TMM Radioactive bound Initial 
(density energy Radius heating cary temperature 
condition 
anomaly) 
1 (REF) 1x108 1(2.5%) 9.21 50 km 0 Reflective 0.72 

2 1x108 — 0.5(1.2596) 9.21 50 km 0 Reflective 0.72 
3 1x108 1.4(3.5%) 9.21 50 km 0 Reflective 0.72 
4 1x108 2(5.0%) 9.21 50 km 0 Reflective 0.72 
5 1x108 1(2.5%) 9.21 25 km 0 Reflective 0.72 
6 1x108 1(2.5%) 9.21 50 km 10 Reflective 0.72 
7 1x108 1(2.5%) 11.51 50 km 0 Reflective 0.72 
8 1x108 1(2.5%) 9.21 50 km 0 Periodic 0.72 
9 1x108 1(2.5%) 9.21 50 km 0 Reflective 0.5 


The radioactive heating rate in case 6 is a non-dimensional value for the TMM materials. The initial temperature is normalized by a reference temperature of 2,500K. 


